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We present a systematic study of the reconstruction of a non-negative function via maximum 
entropy approach utilizing the information contained in a finite number of moments of the function. 
For testing the efficacy of the approach, we reconstruct a set of functions using an iterative entropy 
optimization scheme, and study the convergence profile as the number of moments is increased. We 
consider a wide variety of functions that include a distribution with a sharp discontinuity, a rapidly 
oscillatory function, a distribution with singularities, and finally a distribution with several spikes 
and fine structure. The last example is important in the context of the determination of the natural 
density of the logistic map. The convergence of the method is studied by comparing the moments of 
the approximated functions with the exact ones. Furthermore, by varying the number of moments 
and iterations, we examine to what extent the features of the functions, such as the divergence 
behavior at singular points within the interval, is reproduced. The proximity of the reconstructed 
maximum entropy solution to the exact solution is examined via Kullback-Leibler divergence and 
variation measures for different number of moments. 
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I. INTRODUCTION 

The reconstruction of a non-negative distribution from 
its moments constitutes the so-called classical moment 
problem, and is an archetypal example of an inverse prob- 
lem pQ |2] in mathematical sciences. Owing to its im- 
portance in the context of probability theory and the 
challenging problems of analysis associated with it, the 
moment problem has attracted the attention of a large 
number of researchers from many diverse fields of sci- 
ence and engineering [5lU3|. In the classical Hausdorff 
moment problem (HMP), one addresses the problem of 
reconstructing a non-negative, real valued function f(x) 
in a finite interval [a, b] from a sequence of real num- 
bers. The sequence forms a 'moment sequence' that 
satisfies the Hausdorff conditions [14]. The problem is 
severely ill-posed in the Hadamard sense [15] . For a fi- 
nite number of moments, most of the existing numerical 
methods are susceptible to large instabilities but several 
methods do exist that attempt to construct a regular- 
ized solution by avoiding these instabilities [5] HHi, • The 
HMP has been addressed by using a variety of methods 
(such as Tikhonov's regularization method [T7] and the 
use of Pollaczek polynomial by Viano [HI [19]) but the 
information-theoretic approach is particularly fascinat- 
ing to the physicists. The latter is based on the maxi- 
mum entropy principle (MEP) proposed by Jaynes |20j . 
The MEP provides a suitable framework to reconstruct 
a distribution by maximizing the Shannon information 
entropy [21] and at the same time ensures the matching 
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of the moments of the distribution. Our interest in the 
moment problem stems from the fact that inverse prob- 
lems of this type are frequently encountered in many ar- 
eas of physical, mathematical and biological sciences [5]- 
rHZI E2l |2"5] . A very simple but elegant example is the 
inversion of the specific heat data of solids. It is known 
that the constant volume vibrational specific heat of a 
solid can be expressed as the convolution of a known 
function of the frequency and the vibrational frequency 
distribution function (FDF) [33]. The task of extract- 
ing the FDF by inverting the experimentally-measured 
values of the specific heat at constant volume as a func- 
tion of temperature is a well-known example of an inverse 
problem in solid state physics [23 HI] • 

The focus of our present work is to reconstruct a non- 
negative function very accurately within the framework 
of the MEP from the knowledge of a finite number of 
moments. Although there exists a number of numerical 
procedures that address this problem, most of them be- 
come unreliable when the number of moment constraints 
exceeds a problem-dependent upper limit. A close re- 
view of the methods and the study of the example func- 
tions presented therein immediately reveal the weakness 
of the methods [3 HE]- For example, it is very difficult 
to reproduce accurately the van Hove singularities in the 
frequency distribution of (crystalline) solids or the pres- 
ence of a gap in the density of electronic states in a solid. 
While the algorithm proposed by Silver and Roder [27J 
does reproduce the latter correctly and is capable of deal- 
ing with a large number of moments, we are not aware of 
any systematic study of function reconstruction by this 
approach at this time. It is, therefore, worthwhile to ex- 
plore the possibility of developing a reliable scheme for 
the entropy optimization program and to apply it to a 
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range of non-negative functions having complex struc- 
ture within the interval. 

The rest of the paper is organized as follows. In Section 
II we briefly describe a procedure that has been devel- 
oped recently by us to reduce the moment problem to 
a discretized entropy optimization problem (EOP) [28] . 
We then test our methodology in Section III by examin- 
ing to what extent it is successful in reconstructing a wide 
variety of functions on the basis of input information in 
the form of Chebyshev moments of the functions. The 
convergence behavior of the maximum entropy solution 
is then discussed in Section IV with particular empha- 
sis on the number of moments. The proximity of the 
reconstructed solution to the exact solution for different 
distributions is also studied via Kullback-Leibler [21] di- 
vergence and variation measures [5U] ■ 



II. MAXIMUM ENTROPY APPROACH TO 
THE HAUSDORFF MOMENT PROBLEM 

The classical moment problem for a finite interval [a, 
b], also known as the Hausdorff moment problem [13] , 
can be stated as follows. Consider a set of moments 

Pi = / x l p(x) dx i = 0, 1, 2, . . . , m, i<m (1) 

J a 

of a function p(x) integrable over the interval with fj,j < 
co y x € [a,b]. The problem is to construct the non- 
negative function p{x) from the knowledge of its mo- 
ments. The necessary and sufficient conditions for a solu- 
tion to exist were given by Hausdorff [14 . The moment 
problem and its variants have been studied extensively 
in the literature [T] d [TSl H3 E] • Mead and Papanico- 
laou [9] have, in particular discussed a number of moment 
problems encountered in various of physics. For a finite 
number of moments, the problem is underdetermined and 
it is not possible to construct the unique solution from 
the moment sequence unless further assumptions about 
the function are made. Within the framework of maxi- 
mum entropy principle, one attempts to find a function 
p{x) that maximizes the information entropy functional, 



S[p] 



p{x) ln[p(a;)] dx 



(2) 



subject to the moment constraints defined by Eq. ([I]). 
The resulting solution is an approximate function 
Pme(x), which can be obtained by functional differenti- 
ation of a Lagrangian with respect to the unknown func- 
tion p(x). The Lagrangian is given by, 



L(p,\) = -S[p]+J2*i \ 

i=0 \ Ja 



x n p(x) dx — pi 



(3) 



The normalized function p(x) is often referred to as 
probability density since it is positive semidefinite and 
the interval [a, b] can be mapped onto [0,1] without any 
loss of generality. For a normalized function with po — 
1, the Lagrange multiplier Ao is connected to the others 
via, 



and the maximum entropy (ME) solution can be written 

as, 



Pme(x) = exp \ 



(5) 



where Z is known as the partition function. 

A reliable scheme for handling the entropy optimiza- 
tion problem subject to the matching of the moments was 
discussed by us in Ref. [35J . The essential idea behind the 
approach is to use a discretized form of the Shannon en- 
tropy functional and the moment constraints using an 
accurate quadrature formula. The constraint optimiza- 
tion problem involving the primal variables is then re- 
duced to an unconstrained convex optimization program 
involving the dual variables of the problem. This guaran- 
tees the existence of a unique solution within the frame- 
work of maximum entropy principle. The solution is least 
biased 32] and satisfies the moment constraints defined 
by Eq. ([I]). The procedure consists of: 1) rewriting the 
Lagrangian of the problem in Eq. Q in terms of the dis- 
cretized variables to obtain the ME solution, 2) using the 
resulting ME solution in association with Eq. ([TJ) to re- 
duce the EOP as an unconstrained convex optimization 
problem in dual variables, and finally 3) minimizing the 
objective function in the dual space to obtain the optimal 
solution in the primal space. 

Using a suitable quadrature (e.g. Gaussian) with a set 
of weights Wj's and abscissae Xj's, the discretized La- 
grangian can be written as, 



L(p,X) 



ln 



Pj 




p, 



(6) 

where < p € R n and A = — A E R m , respectively are 
the primal and the dual variables of the EOP. In the 
equation above, we have used the notation pj = uJjPj 
and tij = (xj) 1 . The discretized ME solution is given 
by the functional variation with respect to the unknown 
function as before, 



Now, 



= exp y^tjj K - 1 



j = l,2,...n. 



(7) 



5L 
5p(x) 



= 



Pme{x) = exp l-y^i x l 



(4) Equations (TJ and Q can be combined together and 

the EOP can be reduced to an unconstrained convex op- 
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timization problem involving the dual variables A's: 



mm 



D(X) = 



E 

3=1 



ujj exp 




Ai — 1 — 



^2 Mi 



( 8 ) 

By iteratively obtaining an estimate of A, D(X) can 
be minimized, and the ME solution p(A*) can be con- 
structed from Eq. |7j). The objective function D(X) can 
be minimized by modifying a method, which is largely 
due to Bergman [33] , and was presented and discussed at 
length in Ref. [55] both for the power and the Chebyshev 
moments. For the latter, the ME solution can be shown 
to be expressed in the form of Eq.(|7]) with tij = T*(xj), 
where T*(x) is the shifted Chebyshev polynomials. In 
the following, we apply our algorithm to reconstruct a 
variety of functions corresponding to different number of 
shifted Chebyshev moments. 



III. APPLICATION TO FUNCTION 
RECONSTRUCTION 

We now illustrate the method by reconstructing a num- 
ber of exact functions from a knowledge of their mo- 
ments. For all but one of the examples studied here, 
the moments of the functions can be obtained from an- 
alytical expressions. In the remaining case the moments 
have been calculated numerically using standard double 
precision arithmetic. As mentioned earlier, we map the 
functions onto the interval [0,1] and assume they are nor- 
malized so that the functions can be treated as probabil- 
ity density functions (pdf ) without any loss of generality. 
It is well-known that for a finite number of moments, the 
Hausdorff moment problem cannot be solved uniquely. 
One needs to supply additional information to choose a 
suitable solution from an ensemble of solutions that sat- 
isfy the given moment constraints. The maximum en- 
tropy (ME) ansatz constructs the least biased solution 
that maximizes the entropy associated with the density 
and is consistent with the given moments. The accu- 
racy of the reconstructed solution can be measured by 
varying the number of moments. A comparison with the 
exact solution (if available) would reveal to what extent 
the ME solution matches with the exact solution. For 
an unknown function with a finite set of moments, the 
quality of the ME solution may be judged by the prox- 
imity of the input (exact) moments to the output (ap- 
proximated) moments resulting from the reconstructed 
distribution. By increasing the number of moments one 
can systematically improve the quality of the solution. 
It should, however, be noted, that for a function with a 
complicated structure, the convergence of the first few 
moments does not guarantee its accurate reproduction. 
The ME solution in this case may not represent the ex- 
act solution, but is still correct as far as the maximum 
entropy principle is concerned. It is therefore important 
to study the convergence behavior of the solutions with 



moments for a number of functions with widely differ- 
ent shapes. To this end we compare, in the following, 
our maximum entropy solution corresponding to a vari- 
ety of exact distribution and a distribution amenable to 
an accurate numerical analysis. 



A. Case 1 : f(x) = 1 

We begin with a step function which is unity through- 
out the interval [0, 1]. As mentioned earlier, we use the 
shifted Chebyshev polynomials T*(x), which is defined 



T*(x) = T n (2x-1) 

T n (x) — cos [ncos _1 (a;))] for n = 0, 1, . . . 

The moments can be calculated analytically in this 
case, and are given by, 

1 + (-1)" 

/io = 1; Mi = 0; /i n = - - 2 ^ 2 n,tl. 

Although the function does not have any structure, it is 
particularly important because of its behavior at the end 
points. Owing to the presence of discontinuities at x = 
and 1, the function is difficult to reproduce close to these 
points. The sharp discontinuities cause the reconstructed 
function (from a small number of moments) to exhibit 
spurious oscillations near the end points. The oscillations 
are progressively suppressed by increasing the number of 
Chebyshev moments in our iterative method. Beyond 
100 moments the oscillations completely disappear. This 
behavior is seen clearly in fig [T] where we have plotted the 
reconstructed functions corresponding to 40, 60 and 80 
moments. The oscillations are particularly pronounced 
as one approaches x = 1, but die down with increase 
in the number of moments. The result corresponding 
to 100 moments is presented in fig(5] The plot clearly 
reveals that the function has been reproduced with an 
error, which is less than 1 part in 10 6 . 



B. 



Case 2 : f(x) = 



The next example we consider is a square-root function 
f(x) = | a; 2, where the prefactor is chosen to normal- 
ize the function. In many physical problems, we often 
encounter distributions showing a square-root behavior. 
For example, the spectral distribution of a free electron 
gas in 3-dimension is related to the energy via \f~E^ and 
the square-root behavior persists in the weak interaction 
limit (at low energy). It is therefore important to see 
if such a square-root function can be reproduced with 
a high degree of accuracy using our maximum entropy 
ansatz. The shifted Chebyshev moments for the present 
case are given by, 
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Yin 2 



9 - 40n 2 + 16n 4 



for n > 0. 



The results for the function are plotted in figs. [3] to 
H]for 100 moments. The reconstructed function is found 
to match excellently with the exact function throughout 
the interval as shown in fig(3] Of particular importance 
is the behavior of the function near x — and 1. The 
square-root behavior is accurately reproduced without 
any deviation or oscillation near x = as is evident from 
fig|4j Similarly, the behavior near x = 1 is also repro- 
duced with a high degree of accuracy as shown in fig(5] 
Since our method can exploit the information embedded 
in the higher moments, it is capable of reproducing the 
function very accurately without any oscillation. 



C. Case 3: A double-parabola with a gap 

Having discussed two relatively simple examples, we 
now consider a case where the function vanishes in a fi- 
nite domain within the interval. Such a function appears 
frequently in the context of the energy density of states 
of amorphous and crystalline semiconductors. It is in- 
structive to study whether our maximum entropy (ME) 
method is capable of reproducing a gap in the energy 
eigenvalue spectrum. Since the moments of the electronic 
density of states can be obtained from the Hamiltonian of 
the system, our method can be used as an alternative tool 
to construct the density of states from the moments. This 
is particularly useful for treating a large non-crystalline 
system (e.g. in the amorphous or liquid state), in which 
case the direct diagonalization of the Hamiltonian ma- 
trix is computationally overkill and scales with the cubic 
power of the system size. In contrast, our maximum 
entropy ansatz provides an efficient and accurate proce- 
dure for the determination of total (band) energy and 
the Fermi level subject to the availability of the energy 
moments. Here we use a toy model of a density of states 
that consists of two parabolae separated by a gap to il- 
lustrate the usefulness of our method. In particular, we 
choose a normalized distribution with a gap from x\ to 
x 2 , 



A x (xi — x) 



for x < x\ 



fix) 



B (x — X2) (1 — x) for x > X2, 
where A and B are given by 



A = 



6 

x\{\ +xi — x 2 ) 



D = 



(l-X 2 ) 2 (l+.Tl -X 2 )' 



In the present case, we choose X\ = g and X2 = | giving 
the value of the gap (x 2 — X\) — ^. The Chebyshev mo- 
ments of the function can be calculated exactly, and as in 



the previous examples the function is reconstructed from 
the moments. In figs. [6] and [7] we have plotted the results 
obtained from our ME ansatz along with the exact func- 
tional values at the quadrature points. It is remarkable to 
note that the reconstructed function matches excellently 
with the exact one. Furthermore, the method reproduces 
the gap between the parabolae correctly without any os- 
cillation in the gap. Table 1 lists the size of the gap cor- 
responding to different number of moments for two sets 
of Gaussian points. Since we are using a finite number of 
quadrature points, the accuracy of our gap size is limited 
by the resolution of the (non-uniform) quadrature grid 
near the gap. We have chosen a tolerance e = 5.0 x 10~ 3 
for the reconstructed functional value to locate the onset 
of the gap (i.e. the zero of the function) [34] . It is evident 
from table 1 that as the number of moments increases, 
the size of the gap improves and eventually converges 
very close to the exact numerical value. The accuracy 
can be improved further by using more Gaussian points 
in the quadrature. 



TABLE I: Numerical values of the gap for different number 
of moments from the reconstructed double-parabolic distribu- 
tion. 



Moments 


96 points 


192 points 


20 


0.1622 


0.1676 


40 


0.1813 


0.1875 


60 


0.1821 


0.1902 


80 


0.1823 


0.1909 


100 




0.1922 


Exact numerical 


0.1941 


0.1945 



D. Case 4: f(x) = 



We now consider a function that has singularities in 
the range [0,1]. For the purpose of our discussion we re- 
fer to this function as 'U-function' hereafter. The shifted 
Chebyshev moments of the function have the interesting 
property that except for the zeroth moment, all the other 
moments are identically zero. The task of the ME algo- 
rithm in this case is to construct a function having all 
the moments zero except for the zeroth moment, which 
is unity by normalization. It may be noted that the elec- 
tronic density of states per atom D(E) of an infinite chain 
with a nearest neighbor interaction can be expressed in 
the form, 



D{E) 



1 



1 



7T ^4/3 2 -(E- af 



where a and /3 are the on-site and the nearest neigh- 
bor hopping integrals respectively. The zeroth moment 
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is unity, which implies that there is only one state asso- 
ciated with each atom. For a = and /3 = |, the density 
of states can be mapped onto the U-function within the 
interval [0:1], and our algorithm can be applied to re- 
construct the latter. An important characteristic of the 
density of states (or distribution function) is that it di- 
verges at the band edges (or at the end points). Since 
all the Chebyshev moments are zero aside from the ze- 
roth moment, it is important to see if the algorithm is 
capable of generating the density with the correct diverg- 
ing behavior at the (band) edges. In figj8]we have plot- 
ted the results for the function for three different sets 
of moments M — 10, 40, and 80 to illustrate how the 
approximate solutions improve with the increase of the 
number of moments. The shape of the function begins 
to emerge correctly even for as few as first 10 moments 
but with significant oscillations and poor divergence be- 
havior near the end points. As the number of moments 
increases, the solution rapidly converges and the oscilla- 
tions begin to disappear. In fig|9]we have plotted the re- 
sults for M — 120. The reconstructed function matches 
excellently throughout the interval with the exact one. 
The behavior of f(x) near the left edge at x = is shown 
in fig 10 from £=0 to £=0.05. It is evident from the plot 
that even for very small values of x near the left edge, 
the reconstructed values agree with the exact values ex- 
cellently. A similar behavior has been observed near the 
right edge of the band near x—1. The capability of our 
method in reconstructing a function with singularities in 
the interval is thus convincingly demonstrated. 



E. Case 5: A function with a finite discontinuity 

The functions that we have discussed so far in the ex- 
amples above are continuous within the interval. It would 
be interesting to consider a case where the function has 
a finite discontinuity within the interval. As an example, 
we choose a double-step function, 



/(*) = 



| for x < xi 



§ for x > x\, 



which has a finite discontinuity at x\ = h. It is 
rather challenging to reconstruct the function from the 
moments so that the local behavior near the discontinu- 
ity at x = 1/2 is correctly reproduced. As before, the 
moment integrals can be calculated analytically in this 
case. In fig 11 we have plotted the function for 10, 20 
and 50 moments. The solutions for the first two sets are 
expected to be less accurate, and indeed they show sig- 
nificant oscillations in the figure. For 50 moments the 
match is quite impressive. On adding further moments, 
the solution progressively improves. Figure fl2] shows the 
remarkable accuracy with which the function is repro- 
duced by employing the first 100 moments. An important 



feature of the reconstructed function is that the discon- 
tinuity has been correctly reproduced with the exception 
of two points. From the various cases studied so far, we 
conclude that about 80 to 120 moments are needed for 
point-wise matching of the exact and the reconstructed 
functions. 



F. Case 6 : An unknown density 

Up until now, we have considered cases where the ex- 
act form of the function is known. In practical problems, 
however, it is more likely that the exact function is not 
available. We should therefore consider a case where the 
analytical expression for the distribution is not known, 
but a direct numerical solution is possible. As an exam- 
ple of such a distribution, we choose the natural invari- 
ant density of the logistic map g(x) = V x (1 — x) with 
r = 3.6785. The invariant density for the map can be 
obtained by calculating the moments from the time evo- 
lution of an ensemble of initial iterates Xq as discussed 
in Ref. 35J. Since the map is ergodic for this value of 
r, the moments obtained via the time evolution of the 
map are identical to the moments of the natural invari- 
ant density |23, 35 . The task of our maximum entropy 
algorithm is to reconstruct the approximate density, and 
to compare it with the numerical density. The latter can 
be obtained from a histogram of the iterates and averag- 
ing over a large number of configurations [35] . The result 
from our ME ansatz using the first 80 moments is plot- 
ted in fig j!3| along with the numerical density. The plot 
clearly demonstrates that every aspect of the fine struc- 
ture of the numerical density is reproduced excellently in 
the maximum entropy solution. 

Finally, we end this section considering a rapidly os- 
cillatory function having complex structure within the 
interval [0,1]. An example of such a function can be con- 
structed as, 



f(x) = -(sin(167x) + cos(73x)) + 6 {x - ^) 2 + ^ 



(9) 



where the prefactors are chosen to normalize the func- 
tion. In the context of studying diffusion in a rough 
one-dimensional potential, Zwanzig has studied such a 
function to obtain a general expression for the effective 
diffusion coefficient by analyzing the mean first-passage 
time [36] . The maximum entropy construction of this 



function is plotted in fig 14 for 90 moments along with 
the exact function. Once again the function is reproduced 
excellently with every little details of the local minima 
and maxima of the function. 



IV. CONVERGENCE BEHAVIOR OF THE 
RECONSTRUCTED SOLUTION 

The convergence behavior of the maximum entropy 
solution has been discussed at length in the litera- 
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ture HH [301 EH HB]- The analytical efforts are partic- 
ularly focused on constructing bounds of the proximity 
of the reconstructed density to the exact density assum- 
ing that a given number of moments of the distributions 
are identical. In particular, given a target distribution 
f(x) and a reconstructed distribution /m(^) that have 
identical first M moments (1q = 1, Mi, ■ ■ • , Mm, the prox- 
imity of the two distributions can be expressed via the 
Kullback-Leibler divergence [53] and the variation mea- 
sure l30l. 



D KL [fjM}= f /(x)lnf/M- 
Js 



dx 



dMJm] = I ' \f. 

Js 



m(x) - f(x)\dx, 



(10a) 
(10b) 



where S is the support of the densities f(x) and /m(^)- 
The divergence measure is also known as the relative en- 
tropy or information discrimination, and Dkl > with 
the equality holding if and only if j(x) — juix) for all x. 
A lower bound for the divergence measure Dkl in terms 
of the variation measure was given by Kullback [39] : 



Dkl > 



D 2 



D 4 
12 ' 



(11) 



where it was assumed that the first M moments are iden- 
tical for both the distributions. Since the exact distribu- 
tions are known for the examples considered here (except 
for the case 6), we can use these measures to examine if 
the reconstructed solution indeed satisfies the inequality. 
To this end, we first study the convergence of the mo- 
ments of the reconstructed distributions with iteration 
and establish that the moments can be matched very ac- 
curately so that for practical purposes the reconstructed 
moments can be taken as identical to the exact (input) 
moments for the calculation of measure in Eqs. ( 10a I and 



(10b) 



of the exact moments from the moments provided by the 
ME solution. 



A 1 (N,M) = 



\ 



1 M 



(12) 



Here fi and fi(N, M) respectively denote the exact (or 
input) and the reconstructed (or output) moments, and 
the latter depends on the number of moments (M) and 
the iteration number (N). In the context of our present 
study, the exact moments of the functions are known, 
but in many practical cases they may not be available 
and need to be replaced by the input moments available 
for the problem. The quantity Ai provides a measure of 
the proximity of the first M reconstructed moments to 
the exact ones, and a small value of Ai is indicative of the 
fact that the moment constraints are satisfied with a high 
degree of accuracy. The value of Ai becomes as small as 
10~ 14 provided an adequate number of iterations are per- 
formed to match a given set of moments. In fig |15[ we 
have plotted Ai for the case of U-function with the num- 
ber of iteration progressively increasing to N = 6 x 10 6 . 
The RMS deviation decreases rapidly with the number 
of iteration and eventually drops to a value of the order 
of 10~ 14 . An examination of the data suggests that Ai 
can be fitted to an exponential decay with iteration and 
is plotted in fig 15 This behavior is also observed for the 
other distributions discussed in Section 3. It thus follows 
that the algorithm converges quite rapidly, and that the 
moment constraints can be satisfied to a high degree of 
accuracy even for a very large moment set. 



Convergence with respect to the number of 
moments 



A. Convergence with respect to the number of 
iteration 

As mentioned earlier, we have observed that the qual- 
ity of the ME solutions depend on two factors: 1) the 
number of iterations and 2) the number of moments used 
for the purpose of reconstruction. In general, for a dis- 
tribution with a fine structure, it is difficult to determine 
the minimal number of moments that are needed to re- 
construct the function accurately. However, by studying 
a number of distributions with varying complexities and 
their convergence behavior, it is possible to obtain some 
useful information about the rate of convergence. We ad- 
dress this issue by choosing the U-function (case 4) as an 
example, but the observation is found to be true for other 
cases as well. For a systematic study of convergence be- 
havior of the reconstructed moments with iterations, one 
requires a measure of the goodness of the fit. We there- 
fore introduce Ai, the root mean square (RMS) deviation 



While Ai provides a measure of the goodness of the 
fit for the moments, it does not necessarily guarantee 
a point-wise convergence of the reconstructed function 
with the exact one. This is particularly true if a small 
number of moments are used to reconstruct the function 



that has a fine structure in it (cf. fig 13 ). In this case, the 
reconstructed moments can be matched to a high degree 
of accuracy with input moments, but the solution may 
still miss out the characteristic feature of the distribution 
folded in the higher order moments. The maximum en- 
tropy solution in this case may not reproduce the actual 
solution even though the approximate moments are very 
close to the exact moments. To ensure that Ai indeed 
attains a sufficiently small value, we need to study the 
approximate solution vis-a-vis the number of moments 
for a fixed cycle of iterations. Since the exact functions 
are known in our cases, the simplest way to measure the 
quality of the ME solution is to construct the RMS devi- 
ations of the reconstructed functions from the exact ones 
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in the interval [0,1]: 



A 2 (JV,M) 



\ U 3 i=l 

A 2 (M) for large N, 



(13) 



where n g is the number of points used in the quadrature. 
Here we have assumed that the dependence of A 2 on 
N can be neglected so that fi(N,M) « f%(M), which 
holds for large N owing to the fast decay of Ai. We 
choose Ai = 10~ 15 for each of the moment sets to study 
the variation of A 2 (M) for different values of M. In 
practice, the exact function may not be available but the 
expression for A 2 can still be used by replacing the exact 
function fa by fa (M + AM) and constructing the RMS 
deviation for increasing values M and AM. In fig 16 we 
have plotted A 2 for the case of U-function for different 
values of M. The plot shows a monotonic decrease of A 2 
with the increasing values of M. For this function, we see 
that a value of M=100 to 120 provides a small enough A 2 
to reconstruct the function accurately when Ai = 10 -15 . 
The solid line in the figure is an exponential fit to the 
data indicating a fast convergence of our algorithm with 
respect to the moments for a fixed value of Ai. 

The proximity of the reconstructed distribution to the 
exact one can be quantified in terms of the divergence 
measure and the variation distance as defined in the be- 
ginning of this section. A number of inequalities can be 
found in the literature .'!() 38! that provide lower bounds 
of the relative entropy. The inequality in (11) is an ex- 



ample of such a bound although still sharper bounds are 
available in the literature [30, 38 . In fig 17 we have plot- 



ted the relative entropy of the U-function for different 
number of moments. The reconstructed solution for each 
of the moment sets M corresponds to Ai = 10~ 15 so 
that the first M moments of the exact and the recon- 
structed functions are practically identical to each other. 



The right hand side of the inequality (11 1 is also plotted 



in the same figure for comparison. As the reconstructed 
solution approaches the exact solution with increasing 
number of moments, the relative entropy or information 
discrimination between the two distributions decreases 
and eventually comes very close to the analytically pre- 
dicted lower bound. 



maximum entropy optimization. The method consists of 
mapping the original constraint optimization problem in 
primal space onto an unconstrained convex optimization 
problem in the dual space by using a discrctized form of 
the moments and the Shannon entropy of the function to 
be reconstructed. The resulting optimization problem is 
then solved iteratively by obtaining the optimal set of La- 
grange's parameters as prescribed in Ref . |28) . By virtue 
of its ability to deal with a larger number of moments, 
our present approach is extremely robust and accurate. 
This makes it possible to reconstruct a variety of function 
that are difficult to handle otherwise. 

We demonstrate the accuracy of this method by ap- 
plying to a number of functions for which the ex- 
act moments are available. The method accurately 
reproduces not only smooth and continuous functions 
(such as square-root and double-parabolic functions) but 
also non-smooth and discontinuous functions (such as a 
double-step function with a finite discontinuity). It also 
captures the fine structure in a rapidly oscillatory func- 
tion of known analytical form and the invariant densi- 
ties of a logistic map corresponding to special values of 
the control parameter for which no analytical results are 
available. A convergence study of the reconstructed mo- 
ments suggests that the RMS deviation of the moments 
(from the exact ones) can be made as small as 10~ 15 
indicating the accuracy with which the input moments 
can be matched with the reconstructed ones. A direct 
comparison with the exact functions studied here reveals 
that the method indeed converges to the correct solution 
provided a sufficient number of moments are available as 
input. The general trend of the convergence profile is 
similar in all the cases: the quality of the reconstructed 
function markedly improves with the increase in the num- 
ber of moments. The numerical calculations suggest that 
the convergence toward the exact solution is almost of 
exponential nature. 
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FIG. 1: (Color online) The step function, f(x) = 1, recon- 
structed using the first 40, 60, and 80 Chebyshev moments. 
The plot shows the presence of oscillations at the right edge 
with decaying amplitude as the number of moments increases. 
The data for 40, 60 and 80 moments are indicated in the 
figure by boxes (blue), triangles (red) and circles (green) re- 
spectively. A magnified view of the right edge is shown in the 
inset. 
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FIG. 2: (Color online) The step function as in fig[T] using 
the first 100 Chebyshev moments. The oscillations now com- 
pletely disappear and the function is reproduced with an error 
less than 1 part in 10 6 . 
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FIG. 3: The function, f(x) = %x 5 , and its maximum entropy 
reconstruction using the first 100 Chebyshev moments. For 
clarity, every second data point is plotted in the figure. The 
exact function is evaluated at the quadrature points and is 
drawn as a line for comparison. 
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FIG. 4: The behavior of the function f(x) = for very 
small values of x along with the exact functional values evalu- 
ated at the quadrature points. Note that only one point is off 
the graph indicating an excellent match to the exact function 
for 100 moments. 
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FIG. 5: The behavior of the function f(x) = |cc 2 near x = 1. 
The exact function is also plotted for comparison. 
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FIG. 6: The reconstruction of a function with a gap in the 
interval. The double-parabola with a gap is reconstructed 
using the first 100 moments. The exact values are also plotted 
in the figure for comparison. 
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FIG. 7: The reconstructed double-parabola near the gap 
along with the exact function at the quadrature points. Ow- 
ing to the finite number of quadrature points, the recon- 
structed function has a non-zero value at x% = 0.4 and X2 
= 0.6. 
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FIG. 8: (Color online) Reconstruction of the U-function as 
defined in the text. The data correspond to the first 10, 40 
and 80 moments as indicated in the plot. 
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FIG. 9: The reconstructed U-function using the first 120 mo- 
ments along with the exact function evaluated at the quadra- 
ture points. The reconstructed function matches point- wise 
to the functional values as indicated in the figure. 
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FIG. 10: The divergent behavior of the U-function near x = 
0. The method accurately reproduces the function for very 
small values of x. 
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FIG. 11: (Color online) The reconstructed double-step func- 
tion for the first 10, 20 and 50 moments. The reconstructed 
function improves progressively with the increase in the num- 
ber of moments. 
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FIG. 12: (Color online) The reconstructed double-step func- 
tion using the first 100 moments. The exact function is also 
shown as a line. 
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FIG. 13: A reconstructed density with sharp peaks obtained 
from the first 80 moments. The distribution corresponds to 
the natural invariant density of the logistic map as discussed 
in the text. The line corresponds to the numerical density 
obtained via histogram method. 
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FIG. 14: (Color online) The reconstruction of an oscillatory 
function with a fine structure as discussed in the text. The ex- 
act functional values are also plotted at the quadrature points 
for comparison. The location of the local minima and max- 
ima are excellently reproduced from the first 90 moments of 
the function. 
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FIG. 15: The semi-log plot of the RMS deviation Ai(JV) 
for the U- function with iteration N expressed in unit of 10 6 . 
The RMS values decay exponentially with iteration after an 
initial crossover around N = 0.8. For clarity of presentation, 
every second data point is plotted in the figure. The number 
of moments is indicated in the figure. 
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FIG. 17: The variation of the KL divergence measure (cir- 
cles) of the U-function for different number of moments. The 
function of the variation measure from the right hand side of 
the inequality ( |11[ ) is also plotted for comparison. The data 
correspond to the reconstructed solution with Ai = 1Q _ . 



